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Two-hadron correlations are a sensitive probe of the dynamics of gluon saturation and 
the Color Glass Condensate formalism where the degrees of freedom are Wilson lines of the 
gluon field. It is shown that unlike structure functions in DIS and single hadron production 
in proton-nucleus collisions, higher point functions of Wilson lines appear in two-hadron pro- 
duction cross section. We investigate the energy (x) evolution of these higher point functions 
^Nj ■ using the JIMWLK evolution equations and show that dipole approximation, employed com- 

monly in the literature to fit the di-hadron data measured by the STAR collaboration in the 
forward rapidity region, breaks down. This necessitates an investigation of the full hierarchy 
of the JIMWLK evolution equations for these higher point functions and their solutions. 
This can then be used to make a quantitative analysis of the di-hadron correlations in the 
forward rapidity region in deuteron-gold collisions at RHIC and in the long range rapidity 
correlations observed in proton-proton collisions at the LHC. 
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^ ' 81. Introduction 

^: 

It is an experimental fact, established by HERA experiments on electron-proton 
collisions, that the number of gluons in the wave function of a proton grows fast with 
Bjorken XBj ~ %- where Q 2 is the virtuality of the exchanged photon and S is the 
photon-hadron center of mass energy squared. Bjorken xbj can also be understood 
as the ratio of the energy of the proton carried by a gluon (or any parton in general) 
denoted by x. This came as a surprise and had led to an explosion of ideas on the 
fate of gluons in the limit when x — > 0, also known as the high energy limit of QCD. 

The fast growth of the gluon distribution function can be understood in pertur- 
bative QCD to be due to sequential radiation of large number of soft (in longitudinal 
momentum) gluons which populate the wave function. In other words, the gluon 
splitting function is singular in the limit x — > so that due to the large longitudinal 
phase space available at small x, probability for gluon radiation is large even though 
the coupling constant may still be small. Technically speaking, one needs to re-sum 
quantum corrections which are enhanced by Inl/x which in turn leads to a power 
growth of the gluon distribution function with (l/x) A with A ~ 0.3. 

The power-like rise of the gluon distribution function would eventually lead to a 
power growth of hadronic cross sections which would violate the Froissart bound and 
unitarity. This fast growth however is believed to be tamed by non-linear effects, 
expected to be important when the density of gluons in a proton is so large that the 
probability for gluon recombination becomes of the same order as the bremsstrahlung 
radiation, responsible for the fast growth of the gluon distribution function. 

The gluon recombination effects were first considered by Gribov-Levin-Ryskin 
in a pioneering paper^ where they argued that recombination diagrams should be as 
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important as the bremsstrahlung ones when the gluon distribution function satisfies 
the following condition 

a s xG(x,Q 2 ,b t ) 

SlQ* 1 (M) 

where S± is the unit transverse area. Mueller and QirP then calculated the rele- 
vant diagrams in the double logarithm region and confirmed the GLR expectation. 
Solving equation (|1T|) self-consistently leads to a scale Q 2 s {x) where this relation is 
satisfied. This scale is nowadays called the saturation scale and depends on x as 
well as the impact parameter bf. In case of a nucleus, it will also depend on A, the 
nucleon number, as A 1 / 3 so that one expects that the saturation scale is large at high 
energy (x — > 0) and for large nuclei (A 1 ^ 3 — > oo). This means that this high gluon 
density system is weakly coupled so that one can still use weak coupling techniques 
even though the system is non-perturbative. This is a highly non-trivial observation 
first made by McLerran and Venugopalan. 3 ' 

In this limit the standard expressions for particle production using collinear 
factorization in perturbative QCD breakdown due to two reasons, firstly due to the 
large energy effects which appear as large logs of 1/x and are not re-summed in the 
standard expressions. Secondly, due to the large gluon density the standard twist 
expansion is not valid anymore in the sense that all twist terms are as large as the 
leading twist term. This necessitates introduction of another formalism, capable of 
including both effects. The formalism which generalizes perturbative QCD to include 
both of these effects has come to be known the Color Glass Condensate formalism 
and the high gluon density effects are referred to as gluon saturation. 

§2. The Color Glass Condensate 

In order to go beyond the leading-twist perturbative QCD and to include both 
large logs of energy and gluon density, McLerran and Venugopalan introduced an 
effective action which treats the large x degrees of freedom as color charges p to 
which the gluon fields A 11 couple. Due to the small x approximation, the coupling 
between the color charges p and the gluon field is assumed to be eikonal and given 
by a Wilson line of the gauge field A^. The distribution of color charges p is non- 
perturbative and assumed to be given by a weight functional W[p\. To calculate an 
observable, one solves the classical equations of motion for a fixed color charge p and 
then averages over all color charges p with the weight functional W[p]. The weight 
functions W[p] satisfies the so-called JIMWLK evolution equation^ which describes 
the evolution of W with energy or equivalently with rapidity y (y = log 1/x). The 
JIMWLK evolution equation for rapidity evolution of any operator O can be written 
as 

±V»-k(f (2.1) 
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and U x = U(xt) is the Wilson line in the adjoint representation. This formalism has 
been applied to any high energy process which involves at least one proton or nucleus 
in the initial state. Examples are Deeply Inelastic Scattering (DIS) of electrons 
on protons and nuclei, and particle production in hadronic/nuclear collisions. The 
simplest process to consider is DIS structure functions which are related to the total 
(virtual) photon-hadron cross section. In the rest frame of the target, this process 
can be factored into two parts; first the virtual photon splits into a quark anti-quark 
dipole which then scatters on the target. The probability for the virtual photon to 
split into a quark anti-quark dipole is given by the square of the photon wave function 
and is calculable in QED. The subsequent scattering of the quark anti-quark dipole 
on the target is described by the CGC formalism which describes its energy (rapidity 
or x) dependence. The JIMWLK evolution equation for the rapidity evolution of 
the quark anti-quark dipole can be written as 

& Kf v.) = J ^ (JI,^!,), (* * y- - w« Kf v. tr v, v} 

(2.3) 

where r, s are transverse coordinates of the quark and anti-quark. The left hand side 
of the equation describes the rapidity dependence of the probability for scattering 
of a fundamental (quark anti-quark) dipole on the target proton or nucleus. Due 
to the non-linear term on the right hand side, this probability is unitary in the 
sense that it can never be larger than one, unlike probabilities calculated in leading- 
twist pQCD. There are, however, no known analytic solutions to this equation and 
one has to resort to approximate methods in order to gain further understanding 
of its properties. One very common approximation is the leading N c and dipole 
approximation such that the color average (denoted by the brackets in the above 
equation) of the product of any number of Wilson lines is replaced by products of 
color averages of two Wilson lines. With this approximation, the JIMWLK equation 
for the fundamental dipole (2-pt function) reduces to the Balitsky-Kovchegov (BK) 
equation. 5 ' Normalizing the 2-pt function as 

S(r - s) = i- < tr Vj V s > (2-4) 

the BK equation for the evolution of a fundamental dipole is written as 



d , > N c a s f 2 (r - s 



—S(r - s) = / dTz 

dy 2it 2 J (r — z) 2 (s — z) 



|2 



S(r — s) — S(r — z) S{z — s) 



(2-5) 



The dipole profile function S(r — s) is the building block for many observables such 
as DIS structure functions at small x and forward particle production in proton- 
nucleus collisions at high energy. Quite recently the next-to-leading order corrections 
to the BK equation have been derived and the running coupling solution is obtained 
numerically P The single inclusive hadron production cross section is given by^ 1 ® 



da pA^hX 
dY d 2 P t d?b ~ (2tt) 



1 f 1 dx — \f q/p (x,Q 2 ) S[—P t ,b] D h/q (^,Q 2 ) + 
7r ) Jx F %F y ' X F ,H X 
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f g/p (x,Q 2 ) S A [— P t ,b] D g/h (^-,Q 2 ) 



(2-6) 



where S and Sa are the fundamental and adjoint dipoles. Using phenomenological 
models of the dipole profiles above has led to a quantitative description of the RHIC 
data on forward rapidity hadron production in dA collisions (the result for y = 4 
was a prediction® which was later confirmed experimentally). Most recently the 
solution to the running coupling BK equation has been used to describe the data 
successfully® which gives one more confidence about the applicability of the CGC 
formalism in the forward rapidity region of RHIC. Nevertheless, models based on the 
standard collinear factorization can also describe the forward rapidity data by includ- 
ing shadowing and cold matter energy lossP^ 1 Measurement of other processes, for 
instance photon and dilepton production in the forward rapidity region^ will help 
clarify the dynamics of the forward rapidity particle production (see^ 1 for a review 
of CGC and its applications to particle production in proton- nucleus collisions). 

2.1. Two-particle correlations 

Two-particle correlations are expected to contain more information about the 
dynamics of the process than single inclusive production. In addition to transverse 
momentum dependence of the cross section, one can investigate the angular depen- 
dence of the cross section which can shed more light on the dynamics of the process. 
For instance, in the standard collinear factorization approach to two hadron produc- 
tion, the produced partons are back to back (in Leading Order) so that one expects 
a sharp peak on the away side (n). On the other hand, in the CGC formalism, one 
expects to have a disappearance of the peak, due to shadowing generated by small x 
re-summation and pt broadening due to multiple scattering. The simplest examples 
of two-parton productions in the CGC framework include quark anti-quark^ an( j 
2-gluon 14 ' production in DIS. 

Two-hadron angular correlations in deuteron-gold collisions in the forward ra- 
pidity region have been recently measured by the STAR collaboration at RHIC. For 
central collisions a disappearance of the away side hadron is observed as expected in 
the CGC formalism. The underlying partonic process is a projectile quark scattering 
from the target nucleus and radiating a gluon either before or after the scattering. 
The expressions for this process are given irP^ 1 in momentum space. Later they 
were also derived in the coordinate space which leads to a more compact form of 
the equations.^ We refer the reader td^'El 1 f G r the details of derivation and the 
explicit form of the production cross section. Here we just note that production of 
a quark and gluon in dA collisions involves higher point (more than 2) functions of 
Wilson lines which were not present in DIS structure functions or single inclusive 
hadron production in dA collisions. For instance, the following products of Wilson 
lines, denoted O4 and Oq appear in the production cross section 



Oi(r,f:s) =ti V}t a V f t b [U s 



ab 



1 

N, 



tr Vj V f 



(2-7) 
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O e (r, f : s, s) = tr V r Vj t a t b [U s U ; 



hba 



tr V r V} V-s V} tr V, K f - — tr V r 



where the following identity is used 



(2-8) 



(2-9) 



In principle one will need to know the equation describing the evolution of these 
operators with rapidity, the same way that the JIMWLK-BK equations describes 
the rapidity evolution of the weight function and the 2-pt function. The leading N c 
part of the equation describing the rapidity evolution of 4 Wilson lines is derived irP2J 
while the complete evolution equation for the operators O4 and Oq (which appear in 
the di-jet production) is derived in.^ Here we just write the results and refer the 
reader tcP^ 1 for the details. The evolution equation for the 4-pt function is 



d 2 z( 2 



(r — s) z 



+ 



[r — s) 



(r — z) 2 (s — z) 2 (r — z) 2 {s — z) 2 



Oi(r, r : s) 



1 



(r — s) z 



(r — z) 2 {s — z) 2 
(f — s) 2 



tr Vj V z tr V? Vr tr V} V. + 



^ tr V} V s tr V] V F tr vj V z 

(r — z) z (s — z) z 



(r — s) 



+ 



(r — s) 



r — r) 



(r — z) 2 {s — z) 2 (r — z) 2 {s — z) 2 (r — z) 2 (r — z) 2 
tr Vj V z V? V f V} V s + tr Vj V s V"j V f V? V. 



+ 



1 (r — r) 2 



iv~ (r — z) z {r — zy 



(2-10) 



while the evolution equation for the 6-pt function can be written as 



d_ 

dy 



(0 6 (r,r : s,s)) 



N c a s 
"2(2vr) 2 

(s - s) 2 



d 2 z{ 2 [ 



r — s) 



+ 



r — r 



(r — z) 2 (s — z) 2 (r — z) 2 (r — z) 



+ 



(r — s)" {s — s)" 1 , 1 

( f -z) 2 (s-z) 2 + 3 ( s -z)2(s_ z )2J ^(r, r : s, s) - — 



r — r 



(r — s) z 



(s — rY 



v<+ 7 ^W^ 72 " 7 ^27^ 72] tr V * V ? V ° V s tr V r v z tr V s 

(r — z) z {r — z) z (r — z) z {s — z) z [s — z) z {r — z) ZJ 



+ [ 
+ [ 



(r — z) 2 (f — z) 2 (f — z) 2 (s — z) 2 (r — z) 2 (s — z) 2 ^ T z s s z r s s 



(r — s) 



> + f { \2f S - 2 ,2 ~ ( ( \ 2 ^ 2 , 2 ]tiV r Vr t V s V z UiV z V s UrV s vj 

(r — z) z (s — z) z (s — z) z {s — z) z (r — z) z {s — z) z - > 
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+ [ ( f ~ + il ~ ( r ~ 1 tr V V- ] V U f tr V- U f tr V U- f 

+ 2 7 ^2 tr V r y f V s V s tr V 8 V? tr V z V^ 

(s — z) z (s — z) z 

^j- (f — s) 2 (f — s) 2 (r - s) 2 (r — s) 2 1 



■ (r — z) 2 (s — z) 2 {r — z) 2 {s — z) 2 {r — z) 2 {s — z) 2 (r — z) 2 (s — z) 2 
tr F r U s f tr F f f V s tr U s ^ f 

r (f — s) 2 (r — s) 2 (r — f) 2 i (s — s) 2 -, 

(f — z) 2 (s — z) 2 (r — z) 2 (s — z) 2 (r — z) 2 (r — z) 2 (s — z) 2 (s — z) 2 

tr V r tr V 5 U/ tr V s 

L( r _ Z )2( s _ z) 2 (r _ 2 )2(5 _ z) 2 + (g _ z) 2 { - s _ z) 2\ tV r V r V * V s V z V s V s V z 
+ I(r-Z) 2 (s-Z) 2 (f - Z) 2 (S - Z) 2 ( S - Z )2(s_ z) 2\ tVV r V r VzV s V S V Z V S V S 

r (f — s) 2 (f — s) 2 (r — s) 2 (r — s) 2 

+ L?i~~ 7\27^ _ ^2 ~~ 7i _ Iwi _ -,\2 ~ 7Z~~ 7\2fZ _ ^2 + 



(r — z) 2 (s — z) 2 (r — z) 2 (s — z) 2 (r — z) 2 (s — z) 2 (r — z) 2 {s — z) 2 - 

tr V r Vl V s Vl V z V} V, V} - 2- .„ tr V r V^ V- s V? F s ^ V, ^ 

r (f — s) 2 (f — s) 2 (r — s) 2 (r — s) 2 

^ (f — z) 2 {s — z) 2 (f — z) 2 {s — z) 2 (r — z) 2 (s — z) 2 (r — z) 2 (s — z) 2 - 

trVrV^VsVjVsVj 

_ r (f ~ + (r - ( g - ^) 2 1 tr y r yj 

L (r — z) 2 {s — z) 2 (r — z) 2 {s — z) 2 (s — z) 2 (s — z) 2 

We note that these equations are free of power divergences in the limit the internal 
transverse coordinate approaches any of the external coordinates. Unlike the evo- 
lution equation for the dipole profile (2-pt function) which is rather simple, these 
equations are quite involved and it is not easy to investigate their properties analyt- 
ically (see^J for particular kinematics in which these relations get a bit simplified is 
investigated). In phenomenological applications of CGC to two-hadron production 
it is assumed^ that in the leading N c approximation one can write the product 
of higher point functions as the products of two point functions only (the dipole 
approximation) : 

(0 6 (r,f:M)) - (0 2 {r-s)) (0 2 (r-s)) (0 2 (s-s)) +(0 2 (r-f)) (0 2 (s-s)) (0 2 (s-s)). 

(2-12) 

This assumption has been made by Albacete and MarquetP^ 1 in order to fit the 
data on forward rapidity di-jet correlations in dA collisions at RHIC. Unfortunately, 
this (dipole) assumption is wrong and misses many leading N c contributions to the 
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evolution equation as proven inP^ Therefore, a quantitative understanding of the 
forward rapidity two-hadron correlations from CGC is still lacking. 

There are several issues that need to be understood. First, we note that the 
number (n) of terms in the evolution equations for the higher point functions which 
are formally N c suppressed becomes very large so that when n >> N c , large N c 
approximation may cease to be a good approximation. To see this, it is convenient 
to define normalized operators Sq and S4 such that 

Si(r,f : s) = 1 (O4) 
Ca Cf 

S 6 (r,f:s,s) = -^—{0 6 ) . (2-13) 
La Cf 

Here we focus on Sq since it is the more interesting one. Rewriting the evolution for 
Oq in terms of Sq, we note that the terms in the first few lines in eq. (|2-lip which 
are products of three traces will end up being leading order in N c while the next few 
lines which involve only one trace will be suppressed by N%. We note the number of 
N% suppressed terms is 8 for O4 and ~ 20 for Oq . Making a Gaussian approximation 
(see^ ) to these higher point functions would result in further proliferation of these 
N c suppressed terms. 

The second point that needs to be understood better is the energy dependence 
of expectation value of trace of a large number of Wilson lines that appear on the 
right hand side of eq. (|2Tip . One may expect that these higher point functions will 
grow faster with energy than the two-point function. There is an example where this 
happens, energy dependence of a state of four reggeized gluons has been investigated 
i]J2ni> w here it is found that it has a faster rate of growth with energy than the state 
of two reggeized gluons. Whether a similar thing happens in our case is not known 
for sure but is likely to be true. This would mean that the terms with larger number 
of Wilson lines on the right hand side of evolution equation will grow faster with 
energy than the other terms which have fewer Wilson lines. This stronger energy 
dependence may eventually compensate for the N% suppression. 

Furthermore we note the presence of terms on the right hand side of eq. (|2Tip 
which involve lower point functions, namely the 2-pt and 4-pt functions, reminiscent 
of pomeron loop contribution to the BK equation for the dipole profile.^ We note 
that these terms are in addition to the original 2-pt function which was present in 
the definition of Oq. The origin of these terms seems to be due to kinematics but 
needs to be better understood. Therefore, di-jet production in the forward rapidity 
region of deuteron-gold collisions offer a rich context and a unique opportunity to 
investigate CGC correlations. For this, one needs to solve the full JIMWLK equation 
for the weight function W[p] which can then be used to compute any n-pt function 
in the CGC formalism. 

We note that photon-hadron correlations^ in the forward rapidity region in 
deuteron-gold collisions is another process which is more sensitive to saturation dy- 
namics than single inclusive production. It also has the advantage, compared to 
two-hadron production, that it is sensitive only to the 2-pt function so that one 
could use the latest results for running coupling BK equation in order to make quan- 
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titative predictions. 

Another example of the importance of higher point functions in the CGC frame- 
work is the two-hadron production process in proton-proton or nucleus-nucleus colli- 
sions. In this case one can not solve the problem analytically in the full kinematics of 
CGC and has to resort to numerical methods. Nevertheless, it is common to use the 
so-called kt factorization in the dilute region, i.e., where both the target and projec- 
tile are dilute but the energy of the collision is high enough so that a re-summation 
in x is required (the BFKL region). In this case the two-gluon production cross 
section is given by^* 

I dN 2 \ g 12 f* d 2 h 

\cPpdy p cPqdy q / ~ 64(2vr)6 \haa' Jg'bV Tgcd Jg'dd' ) J \J_ 

\P A{ k 2)P A{ k 4)PA {h)pA (h, 



(p - h) 2 (p - k 2 ) 2 (q - k 3 ) 2 (q - A; 4 ) 2 

x (p* a B ( P ~ k 2 )p* b B (q - k 4 )p B c '(p - k x )p B « (q - k 3 ) ) (2- 1 1 ) 



where p, y p , q,y q are the transverse momenta and rapidities of the two produced 
gluons, A and B label the color charge p of the projectile and target proton or 
nucleus and the Lipatov vertex is denoted by D*. The color charge p is the source 
for the classical gluon field which satisfies the relation (in the covariant gauge) 



A»(x + ,r) = 5 f "-a(x + ,r) = -g 5(x + )^ T p(x + ,r) , (2-15) 



The standard un-integrated gluon distribution function <P(x,p 2 ) is defined in terms 
of the color average of two p's as 

(p* a (k)p b (k>)) (x) = _ M _ T (2vr) 3 5(A ; - k') #(z, k 2 ) (2-16) 



Due to the complexity of the expression, it is common^ to make the dipole approx- 
imation and write (symbolically) 

p a p b p c p d ) = 5 ab 5 cd (p 2 ) 2 + 5 ac 5 bd (p 2 ) 2 + 8 ad 5 bc {p 2 ) 2 + ■■■ (2-17) 



where p 2 denotes (pp). Again, the dipole approximation breaks down under rapidity 
evolution. The correct evolution equation for the product of four p's was derived 
inPS and reads 

(a r a f ot s Ug) = J d z 

(r — z) 2 (f — z) 2 (s — z) 2 (s — z) 2 z 2 
+ l f d 2 z 



it J (2vr) 2 
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(2-18) 



This equation is valid in the dilute region of the proton or nucleus so that all higher 
order terms in the field a have been neglected on the right hand side of the equation. 
Solving this equation (which can be done using numerical methods), one can use 
the solution in the expression for two-gluon production (|2-14|) and investigate the 
dependence of the production cross section on the transverse momenta and rapidities 
of the two gluons. This will be extremely interesting due to the recent results from 
heavy ion collisions at RHIC and proton-proton collisions at LHC where long range 
rapidity correlations are observed!^ These long range rapidity correlations must be 
generated very early after the collision and subsequent re-scattering can not change 
these due to causalityP^ Furthermore, in proton-proton collisions one does not 
expect to create a medium with a size and life time much larger than a Fermi, which 
is how the heavy ion community defines a Quark-Gluon Plasma. Thus, viscous 
corrections to ideal hydrodynamics must necessarily be very large for transverse 
momenta larger than a few hundred MeV. Therefore, one does not expect to see 
flow effects above one 1 GeV . 

In the CGC formalism long range rapidity correlations are due to the production 
of boost-invariant longitudinal color fields at the very early stages of the collision. 
These classical longitudinal fields are boost invariant which leads to production of 
gluons independent of rapidity. Furthermore, the transverse size (correlation length) 
of these fields is inversely proportional to the saturation momentum Q s and de- 
creases with increasing collision energy. This means that most produced particles 
will have transverse momenta of the order of Q s and will be produced approximately 
independently of rapidity. These features of the data from RHIC and LHC are in 
qualitative agreement with the expectations from the CGC formalismp9 Neverthe- 
less, a truely quantitative comparison with the data requires knowledge of the higher 
point functions in CGC formalism. 
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